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1.  Introduction 


Non-crystalline  (amorphous)  ceramics  or  ceramic  glasses  are  used  in  a  variety  of  vital  Army 
personnel,  ground,  and  air  vehicle  applications  that  require  transparent  armor  —  it  is  ubiquitous  in 
tactical  vehicular  windshields  and  side  windows.  Ceramic  glass  is  inexpensive  and  is  formable 
into  large,  flat  plate  and  curved  shapes.  For  many  years  it  has  been  known  that  the  properties  of 
glass  can  be  modified  and  enhanced  through  compositional  modification,  chemical  strengthening, 
annealing,  and  process  control  of  melt  cooling.  Glass  ceramics,  the  controlled  crystallization  of 
nano- sized  single  crystals  in  a  glass  matrix,  offer  another  avenue  for  designed  and  enhanced 
property  modifications  for  transparent  and  opaque  material  applications.  In  addition,  certain 
glass  formulations  have  been  shown  to  exhibit  enhanced  performance  against  shaped-charge  jets 
(SCJs)  (figure  1)  and  other  ballistic  threats,  but  it  is  not  understood  why.  This  is  in  part  due  to 
various  short  and  longer  range  atomic  structural  characteristics  including  atomic  free  volumes, 
cation  coordinations,  bridging  and  non-bridging  oxygen  ( O )  atoms,  bonding  energies,  and 
nanoscale  order  characteristics  (short  and  longer  range)  that  are  difficult  or  impossible  to  quantify 
experimentally  for  ceramic  glass.  In  contrast,  crystalline  ceramics  like  silicon  carbide  ( SiC ), 
aluminum  oxynitride  ( AlON ),  and  others  have  easily  characterizable 

microstructures/mesostructures,  which  consist  of  assemblies  of  individual  single-crystal  grains. 
Ceramic  glasses,  on  the  other  hand,  do  not  have  a  conventional  micro-  or  mesostructure,  as  it  is 
understood  for  crystalline  ceramics.  However,  there  are  micro  structural  scale  variations  in 
ceramic  glass  that  may  include  density  variations  from  atomic  free  volume  variations  or 
microporosity,  size  of  local  atomic  order,  defects  (inclusions,  large  pores,  etc.),  and  others  yet  to 
be  determined.  The  interaction  of  a  stress/shock  wave  from  a  dynamic  impact  involves  many 
structural  changes  not  easily  characterized  by  conventional  equations  of  state  and  can  involve 
reversible  and  irreversible  densification  and  changes  in  bulk  short  range  order  structures 
comparable  to  phase  changes  in  crystalline  ceramics.  For  example,  in  simple  Hertzian 
indentation  testing,  a  wide  range  of  plastic  or  inelastic  deformation  mechanisms  have  been 
observed  in  a  variety  of  glasses.  Multiscale  computational  design  methodologies  (figure  2),  for 
this  class  of  materials  will,  nevertheless,  require  quantitative  and  possibly  statistically  based 
descriptions  of  the  mesoscale,  although  current  efforts  to  develop  such  models  have  fallen  far 
short  of  this  goal. 


1 


(a) 


(b) 


Figure  1.  Enhanced  performance  of  SCJs  into  glass  (a)  test  configuration  for  glass  targets, 
and  (b)  penetration  versus  time  for  several  targets,  after  Moran  et  al.  (i). 
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Figure  2.  A  multiscale  model  for  non-crystalline  ceramics  (glass). 
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Specific  long-term  research  goals  are  threefold: 


1.  Develop  molecular  dynamics  (MD)  process  models  for  a  series  of  chemically  substituted 
amorphous  silica  (a -Si02  or  fused  silica)  materials  for  the  prediction  of  glass  elastic 
properties  assuming  completely  uniform  glass  “mesostructures.”  If  successful,  such  models 
will  enable  ab  initio  prediction  of  structure-property  relations  in  glass  that  will  be  validated 
with  experimentally  determined  elastic  properties. 

2.  Extend  the  MD  models  to  study  densification  of  chemically  substituted  a-,S'/’02  materials 
under  high  pressures  (~60  GPa  Materials  in  Extreme  Dynamic  Environments  [MEDE]) 
relevant  to  ballistic  events  where  reversible  and  irreversible  density  changes  and 
structural  transformations  have  been  observed.  If  successful,  such  models  will  enable  ab 
initio  prediction  of  a -Si02  compressibility,  kinetics,  and  “glass”  phase  transformations 
that  will  be  used  to  develop  equations  of  state  for  a -Si02  materials,  and  thus  form  a 
direct  link  to  the  continuum  scale. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code 
incorporating  the  effects  of  reversible  and  irreversible  densification,  inelastic  deformation, 
and  overlain  with  a  spatiotemporally  evolving  population  of  growing  defects,  which  coalesce 
and  lead  to  ultimate  fracture  and  fragmentation.  It  is  envisioned  that  at  some  time  in  the  not 
too  distant  future,  fully  concurrent  multiscale  computational  finite  element  codes  will  be  used 
by  analysts  on  a  regular  basis  for  optimal  material  design! 

The  remainder  of  the  report  is  organized  as  follows.  General  program  objectives  are  outlined  in 
section  2,  whereas  the  planned  approach  for  modeling  the  multiscale  behavior  of  glass  is 
described  in  more  detail  in  section  3.  A  one-day  short  course  on  “Glass  Technology”  held  at  the 
U.S.  Army  Research  Laboratory  (ARL)  on  October  29,  2010,  is  briefly  described  in  section  4, 
followed  by  an  overview  of  an  ARL/Momentive  Performance  Materials  purchase  order  for 
fabrication  of  pure  and  chemically  substituted  a-S702  materials  that  will  form  the  basis  for  future 
observational  and  computational  efforts  (section  5).  First  principle  quantum  mechanical  methods 
are  used  to  predict  radial  distribution  functions  (RDFs)  in  section  6;  initial  molecular  dynamics 
models  using  simulated  annealing  are  then  used  to  develop  a  glass  shock  Hugoniot  for  fused  silica 
in  section  7.  Initial  efforts  at  computational  modeling  of  plate  impact  experiments  of  fused  silica 
appear  in  section  8.  Section  9  summarizes  the  conclusions  of  the  progress  report. 
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2.  Program  Objectives 


The  long-term  research  goal  of  the  program  is  to  develop  a  concurrent  multiscale  computational 
finite  element  code  for  optimizing  or  enhancing  the  performance  of  various  glasses  against  SCJs; 
the  initial  work  focuses  on  pure  fused-silica  (a -Si02),  and  chemically  varied  a -Si02  materials. 

As  such,  this  objective  falls  squarely  within  the  purview  of  the  Weapons  and  Materials  Research 
Directorate  (WMRD),  since  multiscale  models  are  constitutive  models  (specific  to  a  particular 
material)  wherein  time  evolving  microstructural  changes,  such  as  microcrack  growth,  are  fully 
coupled  to  the  macroscale,  a  phenomenon  that  cannot  be  modeled  or  accounted  for  using  classical 
homogenization  methods.  A  more  immediate  research  objective  is  to  understand  why  certain 
chemically  substituted  a -Si02  materials  exhibit  enhanced  performance  in  the  defeat  of  SCJ  and 
other  ballistic  threats. 

The  program  objectives  are  threefold: 

1.  Develop  MD  process  models  for  a  series  of  chemically  substituted  a -Si02  materials 
for  the  prediction  of  glass  elastic  properties.  This  glass  plays  an  important  role  in  many 
technological  applications  and  its  structure  has  been  inferred  from  neutron-diffraction,  nuclear 
magnetic  resonance,  and  small  angle  x-ray  scattering  (SAXS)  analysis  to  reveal  a  three- 
dimensional  network  consisting  of  tetrahedrally  coordinated  silicon  (Si)  whose  structure  is 
constant  throughout  the  glass  and  defines  its  short  range  order  (SRO).  Long  range  disorder  in 
the  structure  is  manifested  by  a  seemingly  random  variation  in  the  Si-O-Si  bond  angle  in 
adjacent  tetrahedra.  Despite  the  intense  study  of  a -Si02  glass  over  the  last  several 
decades,  much  controversy  still  exists  on  the  best  method  to  model  (i.e.,  via  density  functional 
theory,  molecular  dynamics,  Monte  Carlo  methods,  or  master  equation  techniques)  this 
archetypal  material  for  prediction  of  elastic  properties,  diffusivity,  surface  interactions, 

bond  angle  distribution,  polyamorphism,  and  melt  solidification.  Current  models  that  appear 
in  the  literature  are  often  not  fully  validated  and  progress  towards  this  goal  will  be  made  when 
model  predictions  of  elastic  constants  for  a  series  of  chemically  substituted  a-Si02  glasses 
agree  with  experimentally  determined  constants. 

2.  Extend  the  MD  models  to  study  densification  of  the  chemically  substituted  a -Si02 
materials  under  high  pressures.  Since  long  range  order  in  glass  is  non-existent,  variations  in 
the  SRO,  and  intermediate  range  order  (IRO)  must  be  responsible  for  the  enhanced 
performance  observed  in  ballistic  tests  on  certain  a -Si02  glasses.  If  this  is  the  case, 
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it  may  be  possible  to  use  MD  models  to  predict  macroscopic  ballistic  performance.  Since 
glass  is  subjected  to  extreme  pressure  and  temperature  during  a  SCJ  event,  it  will  be  necessary 
to  study  the  relationship  between  compressibility,  kinetics,  and  phase  transitions  during 
high-pressure  densification  of  a -Si02  glasses  as  manifested  by  changes  in  coordination 
number,  ring  size,  and  free  volume.  Progress  towards  this  goal  will  be  made  when  MD 
derived  equations  of  state  (EOS)  agree  with  those  obtained  experimentally  via  diamond 
anvil  press  and  plate  impact  experiments. 

3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code 
that  incorporates  the  effects  of  reversible  and  irreversible  densification,  and  inelastic 
deformation,  overlain  with  a  spatiotemporally  evolving  population  of  defects  that  grow, 
coalesce,  and  lead  to  ultimate  fragmentation.  This  objective  will  develop  a  computational 
framework  to  combine  the  objectives  from  (1)  and  (2),  and  incorporate  the  influence  of 
fracture  initiation,  growth,  coalescence,  and  fragmentation  of  surface  and  volume  defects 
in  glass  into  a  comprehensive  concurrent  multiscale  finite  element  model  and  code. 
Microcrack  initiation,  growth,  and  coalescence  (sometimes  referred  to  as  failure  waves) 
is  a  multiscale  phenomenon  that  bridges  all  scales  in  a -Si02  glasses  despite  the 
apparent  absence  of  a  structural  mesoscale  for  this  class  of  materials  (see  figure  2). 

Algorithms  for  the  development  of  fully  two-way  coupled  multiscale  codes  are  in  their 
infancy,  and  progress  on  this  objective  will  be  realized  with  successful  development  and 
implementation  of  a  consistent  scheme  for  coarse-graining  localization  phenomena  such  as 
fracture  failure  observed  in  glass. 


3.  Planned  Approach 


The  planned  approach  consists  of  three  components,  which  are  outlined  in  figure  3: 

1.  Validate  the  MD  models  for  a  series  of  chemically  substituted  a -Si02  materials  for  the 
prediction  of  glass  elastic  properties.  Although  there  is  no  effort  within  the  WMRD  to 
predict  a -SiO-2  elastic  properties,  a  hierarchical  multiscale  modeling  effort  is  currently 
underway  within  WMRD,  which  is  focused  on  the  study  of  poly  crystalline  (~200  //m  grain 
size)  AlON  and  validation  of  quantum  and  MD  predictions  of  anisotropic  elastic  constants 
using  diamond  anvil  cell  and  focused-ion-beam  (FIB)/scanning  electron  microscopy  (SEM) 
compression  tests  on  oriented  AlON  single  crystals  (7).  We  plan  to  use  MD  methods 
(with  possible  MD  coarse-graining)  to  simulate  glass  process  modeling  during  melt 
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solidification  by  quenching  a  high-density,  high  temperature,  and  pressure  (with  possible 
polyamorphic  phases)  melt  for  a  series  of  chemically  substituted  a -SiC>2  materials.  Next, 
the  resulting  room-temperature,  chemically  modified  structures  will  be  reversibly  deformed  to 
predict  the  elastic  properties  that  will  be  validated  with  experimentally  determined  elastic 
properties. 


FYll 
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FY14 
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Figure  3.  Five-year  roadmap  consistent  with  the  WMRD  brittle  materials  program. 


2.  Validate  the  MD  models  for  densification  of  chemically  substituted  a -SiC>2  materials 

under  high  pressures.  Although  there  is  currently  no  effort  within  WMRD  to  predict  the  EOS 
of  chemically  substituted  a -S1O2  materials,  MD  methods  have  been  used  to  predict  high- 
pressure  densification  in  these  materials.  MD  simulations  of  pure  a-SiC>2  materials  reveal 
a  Hugoniot  elastic  limit  (HEL)  of  about  10  GPa,  and  an  anomalous  maximum  in 
compressibility  at  around  3  GPa.  Experiments  where  samples  have  been  compressed  to 
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pressures  lower  than  10  GPa  are  indistinguishable  from  the  original  material,  whereas  above 
10  GPa,  materials  can  sustain  an  irreversible  density  increase  from  10-20%  higher  than 
the  starting  material,  although  there  is  controversy  as  to  whether  the  mechanism  is  due  to 
irreversible  coordination  defects  or  permanent  ring  size  modification.  In  contrast  to  the 
behavior  of  a -Si02,  crystalline  quartz  (. a-Si02 ),  undergoes  very  well-known  high 
pressure  polymorphic  phase  transitions  into  a  Coesite  phase  and  a  Stishovite  phase 
(figure  4),  which  involve  changes  in  coordination  of  the  Si  cation  from  four  to  six 
O  atoms.  A  combination  of  diamond  anvil  press  and  plate  impact  experiments 
will  be  conducted  on  a  series  of  chemically  substituted  a -SiC>2  materials  and  compared 
with  MD  densification  simulations  in  glass  in  order  to  understand  the  influence  of  glass 
modifiers  on  changes  in  the  shock  response  of  these  materials.  EOSs  for  a  subset  of 
promising  chemically  substituted  materials  will  be  developed  and  implemented  into  a 
continuum  code  to  determine  if  any  of  the  chemically  substituted  materials  exhibit 
enhanced  ballistic  performance. 


Figure  4.  Crystal  structures  of  quartz,  after  Frye  (2). 


3.  Develop  a  fully  validated  multiscale  finite  element  computational  model  and  code  that 

incorporates  the  effects  of  reversible  and  irreversible  densification,  and  inelastic  deformation, 
overlain  with  an  spatiotemporally  varying  population  of  defects  which  grow,  coalesce,  and 
lead  to  ultimate  material  fragmentation.  The  ultimate  research  objective  is  to  develop  a 
physics-based  multiscale  computational  finite  element  code  for  studying  densification,  and 
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dynamic  fracture  in  non-crystalline  ceramics  (see  figure  2).  Atomistic  behavior 
will  be  linked  to  macroscopic  elastic  properties  and  densification  behavior  through 
development  of  an  EOS  from  first  principles  as  outlined  in  (1)  and  (2)  above.  At  this 
stage,  what  remains  to  be  accomplished,  is  to  successfully  link,  in  a  concurrent  fashion, 
multiscale  failure  phenomena  in  a -Si02  materials  by  incorporating  the  important  role  that 
pre-existing  surface  and  volume  defects  have  on  the  microcrack  growth,  coalescence,  and 
fragmentation  in  this  class  of  materials.  Over  the  past  five  years,  the  first  author  has 
also  been  directly  involved  in  development  of  a  parallel,  concurrent  multiscale  code  for 
heterogeneous  viscoelastic  composites  (8)  under  the  auspices  of  an  ARL/University  of 
Nebraska  cooperative  agreement,  which  will  be  leveraged  and  used  as  the  framework  for  the 
development  of  a  concurrent  multiscale  model  of  a -Si02  materials.  The  chief  challenge  for 
brittle  materials  is  to  correctly  account  for  the  growth  kinetics  of  microcracks  in  a 
multiscale  computational  environment.  The  propagation  of  free  internal  boundaries  at  lower 
scales  will  be  “coarse-grained,”  to  higher  scales  where  global  fracture  failure  and 
fragmentation  is  observed.  As  such,  coarse-graining  algorithms  will  need  to  be  validated 
through  continuum- scale  experiments  on  a-SiO-2  materials  that  measure  dynamic  crack 
propagation  speeds,  mixed-mode  failure,  and  crack  bifurcation  phenomena  using  coherent 
gradient  sensing  and  high  speed  imaging  techniques;  ARL  possesses  capabilities  for 
conducting  such  dynamic  fracture  experiments  in  a -SiC>2  materials  through  the  recent 
establishment  of  a  coherent  gradient  sensing/imaging  facility  funded  by  the  ongoing 
multiscale  modeling  effort  of  AlON .  Models  and  validation  of  the  initiation  and 
propagation  of  discrete  fractures  in  a-SiC>2  materials  should  transition  naturally  into 
models  of  fragmentation  and  comminution  for  behind-armor-debris  applications. 
Fragmentation  experiments  have  classically  been  conducted  using  dynamically  expanding 
ring  experiments  for  defining  fragment  size  versus  strain  rate  and  will  be  used  to 
validate  computational  models  of  fragmentation.  The  development  of  consistent  coarse- 
graining  algorithms  for  fracture  in  materials,  which  is  associated  with  failure  and  loss  of 
material  stability,  is  largely  unexplored  and  is  the  primary  high-risk  goal  of  this  section. 
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4.  Glass  Technology  Short  Course 


A  one-day  short  course  on  the  “Fundamentals  of  Glass  Science,”  was  taught  by  Professor  Arun  K. 
Varshneya,  Alfred  University,  at  ARL  on  October  29,  2010.  Course  attendees  received  copies  of 
Varshneya’s  Fundamentals  of  Inorganic  Glasses  ( 3 ).  The  course  covered  topics  on  (1)  basic 
compositions  and  families  of  commercial  glasses  that  include  vitreous  silicates,  soda-lime 
silicates,  borosilicates,  lead  silicates,  and  aluminosilicates;  (2)  fundamentals  of  the  glassy  state; 
(3)  phase  separation  and  liquid-liquid  immiscibility;  (4)  glass  structures:  oxide  and  non-oxide 
glasses;  (5)  review  of  some  glass  properties  (e.g.,  elastic  properties,  hardness,  viscosity,  thermal 
expansion,  durability);  and  (6)  annealing  and  strengthening. 

The  course  was  well  received,  and  a  question  and  answer  period  ensued  with  questions  like, 
“What  is  the  difference  between  a  glass  and  an  amorphous  material?”  Answer:  “Amorphous 
materials  make  a  continuous  or  discontinuous  volume  transition  to  a  crystal  or  vapor  on  heating, 
whereas,  glasses  continuously  change  to  a  liquid  on  heating.”  See  also,  p.  17  and  18  in 
Varshneya  (3),  and  figure  5. 


Figure  5.  The  volume-temperature  diagram  for  a  glass-forming  liquid,  after 
Varshneya  (3). 
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5.  ARL/Momentive  Performance  Materials,  Inc.  White  Paper  Proposal 


The  following  sections  are  quoted  from  the  Momentive  Performance  Materials,  Inc.*  white  paper 
proposal. 

5.1  Goal: 

The  goal  of  the  work  outlined  in  this  proposal  is  to  produce  several  unique  high  silica 
glasses  and  submit  them  to  the  Army  Research  Laboratory  for  ballistic  performance 
evaluation.  Utilizing  our  internal  capabilities  and  those  of  our  partners  we  will 
propose  select  glass  compositions  and  evaluate  the  glasses  in  terms  of 
thermomechanical,  chemical,  and  mechanical  properties.  A  select  subset  of  glasses 
will  then  be  produced  in  a  form  suitable  for  ballistic  testing.  While  this  brief  study 
will  not  allow  us  to  elucidate  all  the  complex  interactions  between  of  glass 
composition  and  ballistic  mechanical  properties,  this  scoping  trial  should  give  us 
some  guidance  as  to  how  much  additional  improvement  might  be  possible  in 
high-silica,  low  density  glasses  with  additional  research. 

5.2  Background: 

While  there  are  inevitable  tradeoffs  between  weight  and  ballistic  performance, 
improvements  to  the  base  material  can  offer  significant  benefits  to  a  previously 
optimized  package.  For  instance,  fused  silica  has  a  better  stopping  power  and  a 
lower  density  than  soda-lime-silica  glass  resulting  in  a  higher  performance/mass 
figure  of  merit .  However,  fused  silica  suffers  from  a  large  amount  of  secondary 
fracture  created  during  the  initial  impact  making  the  window  less  clear  and  probably 
reducing  it’s  resistance  to  further  impacts.  While  the  cause  of  this  is  not  understood, 
one  might  surmise  that  it  is  driven  by  the  atomic  structure  of  the  glass  (i.e.  number 
of  non-bridging  O,  ring  density,  etc.)  that  can  be  impacted  by  the  chemistry  of  the 
glass.  These  compositional  and  structural  differences  affect  how  the  impact  energy 
transfers  within  the  glass  and  where  the  energy  absorbing  fractures  take  place.  For 
example,  Sehgal  and  Ito  ( 4 )  found  that  changes  in  density  of  silica  glass  had  a  large 
impact  on  brittleness  (see  figure  6  below  from  Sehgal  and  Ito  ( 4 )).  If  a  reduction  in 
brittleness  translates  to  improved  (reduced)  secondary  fracture  and  it  can  be 
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accomplished  without  substantially  increasing  the  density  of  the  glass  (i.e.  staying 
near  high-silica  glass  compositions)  then  significant  optimization  might  be  possible. 


Figure  6.  Brittleness  versus  density  for  glasses  in  the  SiO 2  and  /I^C.j-bascd 
glasses,  after  Sehgal  and  Ito  ( 4 ). 


5.3  Approach: 

We  will  produce  high  silica  glasses  with  additions  of  up  to  10  wt%  of  network 
formers  and  modifiers  using  Momentive’s  lab  scale  puck  fusion  process  at 
temperatures  up  to  1900C.  These  glasses  will  be  characterized  (density,  UV-Vis 
transmission)  and  evaluated  mechanically  to  determine  their  hardness  and 
“brittleness”  following  the  methods  outlined  in  references  ( 4 ,  9),  In  addition, 
viscosity  data  will  be  collected  to  assist  in  selecting  proper  temperatures  for  larger 
scale  melting  of  the  glasses  in  phase  2.  Candidate  compositions  from  phase  one  that 
might  exhibit  improved  secondary  fracture  will  be  melted  in  phase  2  using  our  200  lb 
ingot  test  melter.  These  ingot  samples  will  be  large  enough  to  produce  several 
12"xl2"  windows.  The  cut  and  polished  windows  will  be  furnished  to  ARL  for 
assembly  and  ballistic  testing  to  validate  the  small  scale  mechanical  testing  results 
and  determine  whether  these  relatively  minor  changes  in  composition  can 
significantly  affect  the  secondary  cracking  phenomenon  and  result  in  an  optimized 
high  silica  glass. 

To  date,  fused  silica  samples  have  been  delivered  to  ARL  and  have  been  ballistically  tested. 
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6.  Quantum  Mechanics  Modeling  of  Glass 


We  used  a  first  principles  quantum  mechanics  method,  the  Vienna  ab  initio  simulation  package 
(VASP)  (10)  to  optimize  the  amorphous  structure  of  a  continuous  random  network,  generated  by 
a  Monte  Carlo  bond  switching  model  (11).  A  projector  augmented-wave  basis  with  an  energy 
cutoff  of  600  eV  was  used.  Brillouin-zone  integrations  were  carried  out  using  two  special  k 
points.  Atomic  relaxations  were  performed  until  the  maximum  forces  on  atoms  became  less  than 
0.01  eV/A;  reaction  barriers  were  calculated  using  the  elastic-band  method.  Convergence  tests 
indicate  that  the  energies  reported  have  an  uncertainty  of  less  than  0.05  eV. 

The  optimized  structures  have  different  densities  (2.04-2.25  g/cm3),  different  numbers  of  atoms 
(72,  1 14,  and  192),  and  different  distributions  of  atoms  in  the  rings  (3-7  member  rings  in  the 
structure  shown  in  figure  7(a)  and  4-6  rings  in  the  structures  shown  in  figures  8(a)  and  9(a))  of 
random  continuous  networks.  The  model  has  an  amorphous  structure  inside  the  unit  cell  and 
periodic  boundary  conditions,  and  was  found  quite  efficient  for  probing  of  the  interaction  of 
hydrogen  containing  fractions  with  network  of  silica  (12).  RDFs  (atomic  density  variations  with 
distance  from  a  particular  atom)  calculated  for  these  structures  are  presented  in  figures  7(b),  8(b), 
and  9(b).  Positions  of  two  first  peaks  are  in  reasonable  agreement  with  experimental  RDF 
obtained  from  x-ray  diffraction  (13).  The  calculation  of  bulk  modulus,  B,  from  a  strain-energy 
relationship  and  simple  hydrostatic  pressure  for  a  cubic  unit  cell  with  114  atoms  depicted  in  8(a) 
resulted  in  B  =  40.6  GPa.  We  plan  to  calculate  the  bulk  modulus  and  elastic  constants  for  other 
optimized  models  of  silica. 
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(a) 


(b) 


Figure  7.  (a)  72-atom  unit  cell  with  density  2.15  g/cm3  and  (b)  corresponding  RDF. 


(a)  (b) 

Figure  8.  (a)  1 14-atom  unit  cell  with  density  2.25  g/cm3  and  (b)  corresponding  RDF. 
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(a) 


(b) 


Figure  9.  (a)  192-atom  unit  cell  with  density  2.18  g/cm3  and  (b)  corresponding  RDF. 


7.  Molecular  Dynamics  Modeling  of  Glass 


Time-dependent  processes  for  large  condensed  phase  systems  are  modeled  at  the  atomic  scale 
using  MD,  in  which  classical  equations  of  motion  for  a  system  of  atoms  are  integrated  in  time  (a 
“trajectory”).  Thermodynamic  properties  can  also  be  determined  for  systems  at  equilibrium, 
through  which  time  averages  are  obtained  from  states  sampled  over  the  duration  of  the  trajectory. 
This  method  will  be  used  to  study  the  response  of  glasses  that  have  been  subjected  to  shock 
loading  or  to  and  characterize  properties  that  might  assist  in  understanding  failure  of  these 
materials.  The  MD  method  we  are  using  is  limited  by  the  classical  approximation  and  an 
accurate  description  of  the  interatomic  interactions  for  the  system,  denoted  hereafter  as  the 
potential  energy  surface  (PES).  It  is  particularly  important  that  the  PES  properly  describe  the 
system  under  study  before  it  can  be  used  in  a  predictive  capacity.  Fortunately,  there  are  numerous 
atomistic  models  for  silicate  glasses  that  have  been  developed  for  various  MD  simulation  studies. 
For  this  year’s  effort,  we  used  MD  predictions  of  the  shock  Hugoniot  of  amorphous  silica  for 
verification  and  validation  of  the  various  models,  with  a  goal  to  identify  reasonable  models  for 
use  in  future  MD  explorations  of  material  properties  at  high  loading  rates. 
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7.1  Interatomic  Potentials 


Although  numerous  MD  models  of  silica  and  modified  silicate  glasses  exist,  we  have  focused  on 
two  interatomic  potentials  for  silica.  These  were  selected  for  initial  exploration  because  of  their 
extensive  validation  over  a  wide  range  of  thermodynamic  states  and  because  they  had  parameters 
for  the  chemically  substituted  systems  of  interest  (i.e.,  calcium  [Co],  sodium  [N o] ,  and  aluminum 
[ Al ]).  The  first  potential  has  the  following  form 


u (r)  =  +  Dij  [{1  -  exp  {-al3  •  (r  -  r0))}2  -  l]  +  ,  (1) 

where  k,  qi,  qj,  DtJ,  atj,  r0,  and  C,3  are  constants,  which  are  provided  in  Pedone  et  al.  (14).  The 
first  term  is  a  long-range  Coulomb  interaction,  which  is  accurately  evaluated  through  the  Ewald 
summation  (15)  in  MD  simulations  assuming  a  finite  simulation  cell  with  periodic  boundary 
conditions  imposed,  the  second  term  is  a  short-range  Morse  function,  and  the  last  is  a  purely 
repulsive  term.  For  simplicity,  we  will  refer  to  equation  1  as  a  Morse  potential.  A  second 
potential  was  also  used  in  this  study,  a  Buckingham  potential  in  the  form 

m  =  ® 

where  k,  qi,  qj,  Al3,  rl3,  p,  and  CtJ  are  constants,  which  are  provided  by  van  Beest  et  al.  (7<5);  we 
refer  to  equation  2  as  the  BKS  potential. 

7.2  Generating  Equilibrated  Initial  Structure  of  Fused  Silica 

We  started  with  1536  molecules  (1536  Si  atoms  and  3072  0  atoms)  of  quartz  crystal  in  the 
/3-cristobalite  phase.  This  is  a  tetragonal  phase;  in  order  to  generate  a  cubic  simulation  cell,  the 
z-dimension  was  reduced  to  equal  the  x-  and  y-dimensions,  and  the  dimensions  of  the  simulation 
cell  are  chosen  so  that  the  density  of  the  material  is  equal  to  the  experimental  value  of  2.20  g/cm3. 
This  compression  introduces  considerable  stress  to  the  system;  however,  this  is  not  a  concern  as 
very  high  stresses  are  obtained  during  the  annealing  process. 

To  generate  an  amorphous  phase,  we  anneal  a  system  from  a  melt  by  decreasing  the  temperature 
in  discrete  intervals.  The  system  is  simulated  in  the  constant  volume-constant  temperature 
(NVT)  ensemble  to  ensure  the  density  remains  constant.  We  used  the  molecular  dynamics  code 
LAMMPS  (17),  which  is  equipped  to  support  this  annealing  method  within  a  single  simulation 
run. 
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Using  the  Morse  potential  (equation  1),  the  temperature  was  decreased  in  discrete  intervals 
following  the  schedule  provided  in  Pedone  et  al.  (14).  Starting  at  5000  K,  the  system  was 
equilibrated  by  scaling  the  velocity  every  time  step  for  6000  steps,  then  scaling  was  reduced  to 
every  40  time  steps  for  6000  steps,  and  finally  velocity  scaling  was  turned  off  for  the  final  8000 
steps.  The  temperature  of  the  system  was  then  decreased  by  500  K,  and  the  process  was  repeated 
until  the  temperature  of  the  system  reached  300  K.  This  results  in  a  total  simulation  time  of  470 
ps,  which  is  a  nominal  cooling  rate  of  10  K/ps. 

7.3  Characterizing  Fused  Silica 

The  final  configuration  from  the  annealing  simulation  described  in  the  previous  section  was  used 
as  the  initial  configuration  in  constant  pressure-constant  temperature  ( NPT )  simulations  using 
both  forms  of  the  interatomic  potential  (equations  1  and  2).  All  simulations  in  this  section  were 
run  using  the  DL_POLY  MD  suite  of  codes  (18).  The  timestep  for  the  NPT  simulations  was  0.2 
fs,  with  barostatting  and  thermostatting  rates  of  0.5  and  1.0  ps,  respectively.  The  simulations 
were  run  for  60  ps,  following  a  10-ps  equilibration  with  velocity  scaling  every  10  time  steps.  The 
results  were  analyzed  to  obtain  structural  details  of  the  amorphous  material. 

Our  first  criterion  in  evaluating  the  interatomic  potential  is  the  ambient  state  density,  which  we 
compare  to  the  experimental  value  of  2.20  g/cm3.  The  Morse  potential  resulted  in  a  density  of 
2.26  g/cm3,  slightly  higher  than  the  experimental  value,  while  the  BKS  potential  produced  an 
even  higher  density  of  2.34  g/cm3.  The  density  of  the  Morse  potential  is  within  3%  of  the 
experimental  value,  but  the  error  in  the  BKS  density  is  greater  than  6%  ,  which  is  higher  than 
typically  acceptable.  We  also  calculated  partial  RDFs  for  the  three  atomic  pairs,  Si-O,  Si-Si  and 
O-O,  using  both  potentials.  The  difference  in  RDFs  between  the  Morse  and  BKS  potentials  was 
negligible;  hence,  in  figure  10  we  compare  the  Morse  potential  RDFs  to  those  from  a  shell  model 
potential  for  silica  (5).  Except  for  a  slight  difference  in  the  height  of  the  first  peak,  figures  10(a) 
and  10(b)  are  nearly  identical.  This  indicates  that  the  Morse  and  BKS  models  accurately  capture 
the  expected  structure  of  fused  silica.  We  note  that  the  Tilocca  et  al.  (5)  model  was  not 
considered  for  use  in  this  project  due  to  numerical  difficulties  associated  with  integrating 
trajectories  using  this  model  under  compression. 

Finally,  we  calculate  the  Hugoniot  curve  for  the  two  potentials,  following  the  method  of 
Erpenbeck  (19).  The  Hugoniot  function  is 

H(T,  V)  =  E(T ,  V)  -  E0(T0,  V0)  +  ±  [. P(T ,  V)  +  P0(T,  V )]  (V  -  V0)  ,  (3) 
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Figure  10.  RDFs  of  (a)  our  system  using  the  Morse  potential  and  (b)  system  in 
Tilocca  et  al.  (5). 


where  T  and  P  represent  the  temperature  and  pressure,  V  and  E  represent  the  volume/unit  mass 
and  internal  energy/unit  mass,  and  the  subscript  0  indicates  the  unshocked  reference  state.  A 
Hugoniot  curve  represents  a  set  of  (P,  V,  T )  points  for  which  equation  3  is  equal  to  zero.  To 
calculate  a  single  point,  a  series  of  constant  volume-constant  temperature  (NVT)  simulations  is 
performed  over  a  range  of  temperatures,  all  at  the  same  volume.  The  Hugoniot  as  a  function  of 
temperature  is  fitted  to  a  polynomial,  and  the  H  =  0  value  is  interpolated.  (For  our  calculations, 
a  linear  fit  was  sufficient  to  extract  the  Hugoniot  temperature).  Pressure  data  as  a  function  of 
temperature  can  then  be  fitted  to  a  polynomial  to  determine  the  corresponding  pressure  at  the 
Hugoniot  temperature.  Figure  1 1  shows  the  Hugoniot  curves  for  both  the  Morse  and  BKS 
potentials,  represented  as  pressure  versus  volume  ratio  (V /Vq  or  relative  volume),  as  well  as  the 
experimental  values  for  comparison.  We  see  that  at  high  values  of  V /Vq,  the  Morse  potential  is 
nearly  in  agreement  with  experimental  results,  while  at  low  values  of  V /Vq  the  pressure  is  too 
high.  Conversely,  the  BKS  potential  shows  better  agreement  at  low  values  of  V/Vq,  but  the 
pressure  is  too  low  at  higher  values  of  V /Vq.  We  performed  several  numerical  experiments  in 
which  we  scaled  various  parameters  of  the  two  PES  functions  in  order  to  obtain  good  agreement 
of  the  shock  Hugoniot  at  both  high  and  low  densities;  however,  any  improvement  in  the  Hugoniot 
curve  corresponded  to  an  increased  error  in  the  density.  The  discrepancies  between  our  Hugoniot 
curves  and  the  experimental  values  render  these  models  inadequate  for  our  modeling 
requirements.  We  are  continuing  to  evaluate  other  a -SiC>2  atomistic  models,  as  well  as  exploring 
the  possibility  of  fitting  parameters  using  force  matching  techniques  of  Izvekov  et  al.  (20). 
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Figure  11.  Comparison  of  Hugoniot  curve  of  fused  silica 

(a -Si02)'-  experiment  (•),  Morse  potential  (x)  and 
BKS  potential  (A). 


8.  Equation-of-state  Modeling  of  Glass 


Validation  of  the  equation  of  state  (EOS)  constructed  from  multiscale  modeling  will  be  through 
comparison  with  continuum  data.  Bulk  modulus  data  exist  in  the  literature  for  several  common 
S'f02-based  glasses  measured  in  high  pressure  fixtures  and  diamond  anvil  cells  (DACs)  for 
pressures  up  to  approximately  10  GPa.  Similar  quasi- static  data  can  be  extracted  for  the 
compositions  examined  in  this  study  if  resources  permit.  In  the  neighborhood  of  approximately 
10  GPa  pressure,  many  glasses  begin  permanent  rearrangement  on  the  atomic  scale.  While  it  is 
possible  to  continue  measurement  of  bulk  modulus  in  a  DAC  beyond  this  point,  there  is  no 
established  way  to  measure  density  in  the  amorphous  glasses.  Consequently,  validating  the 
model  predictions  at  higher  pressure  with  quasi- static  data  will  be  challenging.  Experimental 
techniques  development  is  needed  to  pursue  validation  at  higher  quasi-static  pressures. 
Alternatively,  EOS  model  validation  at  high  pressures  can  also  be  accomplished  through 
comparison  with  Hugoniot  data  from  shock  experiments.  One  must  also  be  careful  when  using 
the  shock  data,  but  the  approximations  in  the  data  interpretation  tend  to  diminish  as  pressure 
increases.  Two  of  the  features  that  must  be  dealt  with  are  strength  and  non-recoverable  volume 
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change.  To  determine  the  extent  to  which  these  features  complicate  the  data,  finite  element 
simulations  replicating  symmetric  impact,  gas  gun  Hugoniot  experiments  were  run  in 
ALE3D  (27)  using  the  Livermore  Equation  of  State  (LEOS)  tables  for  fused  silica.  The 
simulations  were  run  with  and  without  strength  to  quantify  the  differences.  The  constant  strength 
model  is  used  to  introduce  a  limit  to  the  elastic  shear  response.  Data  on  inelastic  shear  strength 
as  a  function  of  deformation  were  not  found  in  the  literature,  so  the  constant  flow  strength  is  a 
crude  approximation  for  a  behavior  that  likely  depends  on  density,  strain,  strain  rate,  pressure,  and 
temperature.  Figure  12  presents  the  axial  stress  distribution  at  the  shock  front  for  an  impact 
velocity  of  4  km/s.  Shown  is  the  expected  three-wave  structure  (traveling  to  the  right  in  the  plot), 
with  the  first  step  due  to  the  limit  on  the  elastic  shear  stiffness  and  the  second  due  to  the  volume 
collapse. 
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Figure  12.  Stress  wave  distribution  in  fused  silica. 


The  wave  front  is  not  steady  at  this  lower  impact  pressure  in  that  the  three-wave  fronts  continue  to 
separate  over  time.  A  steady  wave  is  assumed  when  analyzing  data  in  the  context  of  the 
Hugoniot  equations,  and  that  introduces  uncertainty  when  computing  the  energy  dissipation  and 
stress.  One  is  faced  with  the  question  of  which  wave  front  to  use  to  compute  the  wave  speed 
when  analyzing  the  data.  The  experimental  results  under  these  low  pressure  conditions  must  be 
put  in  the  context  of  how  the  data  were  reduced.  At  higher  impact  pressures,  the  plastic  wave  and 
the  transformation  front  both  overtake  the  elastic  precursor  and  create  a  strong  shock.  Under 
these  higher  pressure  conditions,  the  data  analysis  is  not  ambiguous. 
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The  plot  in  figure  13  provides  some  indication  of  the  error  in  the  energy  calculation.  The  error 
metric  plotted  is  the  difference  between  the  energy  computed  by  the  finite  element  code  and  the 
ideal  energy  computed  assuming  the  Hugoniot  relations  apply.  Normalization  is  by  the  ideal 
energy.  It  can  be  seen  that  the  deviation  from  the  Hugoniot  conditions  is  significant  at  lower 
compressions  but  nearly  zero  at  densities  above  4.3  g/cm3.  This  coincides  with  development  of  a 
strong  shock  as  described  previously.  It  should  also  be  noted  that  the  energy  error  is  significantly 
greater  for  the  calculation  with  strength.  The  energy  from  inelastic  dissipation  is  not  captured 
properly  by  the  Hugoniot  energy  balance  in  the  low  pressure  regime  where  the  elastic  precursor 
outpaces  the  plastic  wave.  The  error  vanishes  when  the  plastic  wave  overtakes  the  elastic 
precursor  and  the  wave  becomes  steady. 
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Figure  13.  Energy  error  vs.  density  in  fused  silica. 
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The  stress  density  plots  from  the  simulations  are  shown  along  with  data  by  Marsh  (6)  in  figure 
14.  It  can  be  seen  that  the  LEOS  table  reproduces  the  data  reasonably  well  in  the  high  pressure 
range.  The  strength  used  in  the  simulations  can  make  an  appreciable  difference  at  lower  pressure 
but  because  the  strength  was  specified  as  constant  rather  than  pressure  dependent,  the  relative 
influence  decays  as  the  pressure  increases.  Simulations  were  not  run  at  lower  pressures  because 
the  results  are  increasingly  dominated  by  the  assumed  strength  model  and  volume  collapse, 
neither  of  which  is  thought  to  be  accurate  in  these  simulations.  The  initial  density  is  2.204  g/cm3, 
so  there  is  a  significant  span  of  density  change  where  the  volume  collapse  and  strength  dominate. 
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These  simple  analyses  illustrate  that  strength  effects  on  gas  gun  results  for  glass  can  be 
appreciable,  because  the  strength  is  a  significant  fraction  of  the  total  stress.  The  strength  and 
whether  or  not  the  waves  are  steady  should  be  considered  when  using  the  data  to  calibrate  an  EOS 
or  validate  multiscale  results. 
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Figure  14.  Stress  wave  vs.  density  in  fused  silica 
compared  with  data  from  Marsh  (6). 


9.  Conclusions 


This  six-month  progress  report  on  multiscale  modeling  of  noncrystalline  ceramics  (glass)  has 
focused  on  establishing  the  framework  for  development  of  a  multiscale  computational 
methodology  for  optimizing  or  enhancing  the  performance  of  fused  silica  materials  not  yet 
synthesized.  A  more  immediate  research  objective  is  to  understand  why  certain  chemically 
substituted  fused  silica  materials  exhibit  enhanced  performance  in  the  defeat  of  SCJs  and  other 
ballistic  threats.  Conclusions  consistent  with  the  milestones  shown  in  the  five-year  roadmap 
shown  in  figure  3  are  as  follows: 

1.  Professor  Arun  K.  Varshneya  of  Alfred  University,  gave  a  short  course  on  “Fundamentals  of 
Glass  Science,”  at  the  ARL  on  October  29,  2010;  course  attendees  received  copies  of 
Fundamentals  of  Inorganic  Glasses  ( 3 ),  as  well  as  a  variety  of  practical  considerations 
regarding  glass  formation  and  technology. 
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2.  Fused  silica  specimens  have  been  delivered  to  ARL  under  the  auspices  of  an  ongoing 
cooperative  agreement;  these  specimens  have  been  ballistically  tested,  and  will  serve  to 
validate  future  computational  models  of  fused  silica. 

3.  A  first  principles  quantum  mechanics  method  (VASP)  (10)  was  used  to  optimize  the 
amorphous  structure  of  a  continuous  random  network,  generated  by  a  Monte  Carlo  bond 
switching  model  (11).  A  projector  augmented-wave  basis  with  an  energy  cutoff  of  600  eV 
was  used.  The  optimized  structures  have  different  densities  (2.04-2.25  g/cm3),  different 
numbers  of  atoms  (72,  114,  and  192),  and  different  distributions  of  atoms  in  the  rings  (3-7 
member  rings  in  the  structure  shown  in  figure  7(a)  and  4-6  rings  in  the  structures  shown 

in  figures  8(a)  and  9(a))  of  random  continuous  networks.  RDFs  calculated  for  these 
structures  are  presented  in  figures  7(b),  8(b),  and  9(b).  Positions  of  two  first  peaks 
are  in  reasonable  agreement  with  experimental  RDF  obtained  from  x-ray  diffraction  (13). 

4.  Initial  fused  silica  glass  structures  were  generated  via  LAMMPS  (17)  simulations  by 
simulated  annealing  the  system  from  a  melt  via  a  Morse  potential  (equation  1)  using  the 
schedule  provided  in  Pedone  et  al.  (14).  Final  configurations  from  the  simulated  annealing 
runs  were  used  as  the  initial  configurations  for  DL_POLY  (18)  predictions  of  RDFs  for  fused 
silica,  which  compared  favorably  with  the  results  of  Tilocca  et  al.  (5); 

5.  Hugoniot  curves  for  fused  silica  were  calculated  using  the  method  of  Erpenbeck  (19)  using 
both  Morse  and  BKS  potentials.  The  Morse  potential  accurately  models  pressure  at  a  relative 
volume  greater  than  0.7,  yet  overpredicts  pressure  at  a  lower  relative  volume.  The  BKS 
potential  accurately  predicts  pressure  at  a  relative  volume  on  the  order  of  0.5,  but 
underpredicts  pressure  at  larger  relative  volumes.  The  discrepancies  between  the  computed 
Hugoniot  values  and  experimental  data  render  these  potentials  inadequate  for  our  modeling 
requirements,  and  we  are  evaluating  other  a -Si02  atomistic  models,  as  well  as  exploring 

the  possibility  of  fitting  parameters  using  force  matching  techniques  of  Izvekov  et  al.  (20). 

6.  Symmetric  gas  gun  impact  simulations  on  fused  silica  were  conducted  using  the  ALE3D  code 
based  on  data  from  the  LEOS  tables.  The  classic  three-wave  structure  was  observed  in  the 
fused  silica  simulations,  and  stress  versus  density  calculations  reproduced  the  experimental 
data  reasonably  well  in  the  high  pressure  range.  The  strength  used  in  the  simulations  can 
make  an  appreciable  difference  at  lower  pressure  because  the  strength  was  specified  as 
constant  rather  than  pressure  dependent,  and  it  was  shown  that  the  relative  influence  decays 
as  the  pressure  is  increased. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


a-SiO  2 

crystalline  quartz 

Al 

aluminum 

ALE3D 

arbitrary  Lagrangian-Eulerian  three  dimensional  hydrocode 

AlON 

aluminum  oxynitride 

ARL 

U.S.  Army  Research  Laboratory 

a -Si02 

fused  silica  or  amorphous  quartz 

BKS 

MD  potential  named  after  authors  (16) 

Ca 

calcium 

DAC 

diamond  anvil  cell 

DL_POLY 

MD  code  named  after  developers  at  Daresbury  Laboratory 

DSI 

Director’s  Strategic  Inititiative 

EOS 

equation  of  state 

FIB 

focused-ion-beam 

HEL 

Hugoniot  elastic  limit 

IRO 

intermediate  range  order 

LAMMPS 

large-scale  atomic/molecular  massively  parallel  simulator 

LEOS 

Livermore  equation  of  state 

MD 

molecular  dynamics 

Na 

sodium 

MEDE 

materials  in  extreme  dynamic  environments 

NPT 

constant  pressure-constant  temperature 

NVT 

constant  volume-constant  temperature 

0 

oxygen 

PES 

potential  energy  surface 

RDF 

radial  distribution  function 

SAXS 

small  angle  x-ray  scattering 

SCJ 

shaped  charge  jet 

SEM 

scanning  electron  microscopy 

Si 

silicon 

SiC 

silicon  carbide 

SRO 

short  range  order 

VASP 

Vienna  ab  initio  simulation  package 

WMRD 

Weapons  and  Materials  Research  Directorate 
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RDRL  WMM  B 
T  BOGETTI 
B  CHEESEMAN 
C  FOUNTZOULAS 
A  GAZONAS 
G  GAZONAS 
D  HOPKINS 
PMOY 


B  POWERS 
C  RANDOW 
TSANO 

M  VANLANDINGHAM 
R  WILDMAN 
CF  YEN 
RDRL  WMM  C 
J  LA  SCALA 
RDRL  WMM  D 
ECHIN 
KCHO 

RDRL  WMM  E 
J  ADAMS 
M  COLE 
TIES  SEN 
J  LASALVIA 
P  PATEL 
J  SANDS 
RDRL  WMM  F 
L  KECSKES 
H  MAUPIN 
RDRL  WML  G 
J  ANDZELM 
A  RAWLETT 
RDRL  WMP 
S  SCHOENFELD 
RDRL  WMP  B 
R  BECKER 
S  BILYK 
D  CASEM 
J  CLAYTON 
M  GREENFIELD 
C  HOPPEL 
R  KRAFT 
B  LEAVY 
BLOVE 

M  RAFTENBERG 
S  SATAPATHY 
M  SCHEIDLER 
T  WEERASOOYIA 
RDRL  WMP  C 
T  BJERKE 
S  SEGLETES 
RDRL  WMP  D 
R  DONEY 
D  KLEPONIS 
J  RUNYEON 
B  SCOTT 
H  MEYER 
RDRL  WMP  E 
M  BURKINS 
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NO.  OF 
COPIES 


NO.  OF 

ORGANIZATION  COPIES  ORGANIZATION 

RDRL  WMP  F 
M  CHOWDHURY 
A  FRYDMAN 
N  GNIAZDOWSKI 
R  GUPTA 

RDRL  WMP  G 
N  ELDREDGE 
D  KOOKER 
S  KUKUCK 
G  R  PEHRSON 
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